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Abstract.  In  this  study,  we  describe  the  algebraic  computations  required  to  implement  the 
stochastic  finite  element  method  for  solving  problems  in  which  uncertainty  is  restricted  to  right 
hand  side  data  coming  from  forcing  functions  or  boundary  conditions.  We  show  that  the  solution 
can  be  represented  in  a  compact  outer  product  form  which  leads  to  efficiencies  in  both  work  and 
storage,  and  we  demonstrate  that  block  iterative  methods  for  algebraic  systems  with  multiple  right 
hand  sides  can  be  used  to  advantage  to  compute  this  solution.  We  also  show  how  to  generate  a 
variety  of  statistical  quantities  from  the  computed  solution.  Finally,  we  examine  the  behavior  of 
these  statistical  quantities  in  one  setting  derived  from  a  model  of  acoustic  scattering. 
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1.  Introduction.  It  is  common  practice  for  mathematical  models  to  be  studied 
under  the  assumption  that  data  defining  the  models  are  precisely  understood.  In  real¬ 
ity,  however,  this  simplifying  assumption  is  often  not  valid,  and  there  is  considerable 
uncertainty  in  specification  of  models.  Sources  of  uncertainty  include  geological  prop¬ 
erties  of  transporting  media,  material  properties  of  structures,  and  unknown  aspects 
of  boundary  conditions. 

One  approach  for  addressing  this  issue  is  to  treat  poorly  specified  data  as  random 
variables  having  some  given  statistical  properties  such  as  means  and  higher  order 
moments,  and  then  to  determine  analogous  statistical  properties  of  solutions.  For 
boundary  value  problems  with  uncertain  data  ( stochastic  partial  differential  equa¬ 
tions ),  a  methodology  known  as  the  stochastic  finite  element  method  has  generated 
considerable  attention  in  the  last  decade  [6,  7,  11,  12,  15].  This  approach  starts  with  a 
boundary  value  problem  in  d-dimensional  physical  space.  The  stochastic  component 
of  the  problem  statement  is  then  specified  approximately  using  an  m-dimensional  aux¬ 
iliary  space  which  is  derived  from  an  underlying  probability  space  associated  with  the 
data.  The  result  is  a  (d  +  m)-dimensional  model,  which  can  be  stated  in  a  weak  form 
on  a  suitable  function  space  using  a  combination  of  standard  variational  constructions 
for  the  physical  component  of  the  problem  together  with  averaging  for  the  stochastic 
component.  We  will  outline  the  details  of  this  methodology  in  Section  2. 

Once  this  weak  formulation  is  specified,  a  numerical  solution  of  the  stochastic 
partial  differential  equation  can  be  computed  in  essentially  the  same  manner  as  for 
deterministic  problems.  In  particular,  the  introduction  of  finite  dimensional  subspaces 
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leads  to  an  algebraic  system  of  equations  whose  solution  can  be  used  to  approximate 
statistical  properties  of  the  physical  solution,  such  as  its  mean,  variance  and  covari¬ 
ances.  Our  concern  in  this  paper  is  to  explore  the  computational  costs  of  solving  the 
systems  in  question  and  of  generating  statistical  analyses  of  the  solution. 

We  will  focus  on  problems  where  randomness  only  affects  the  right  hand  sides  of 
the  algebraic  systems,  that  is,  where  the  forcing  terms  or  boundary  data  are  random 
functions.  A  natural  example  of  this  arises  in  models  of  acoustic  or  electromagnetic 
scattering,  where  lack  of  information  about  the  material  properties  of  scatterers  or 
the  shape  and  structure  of  boundaries  such  as  ocean  bottoms  leads  to  uncertainty 
in  boundary  conditions.  We  will  use  this  model,  specifically,  the  numerical  solution 
of  the  Helmholtz  equation,  as  a  benchmark  problem,  and  in  our  assessment  we  will 
explore  computational  issues  associated  with  quantities  such  as  moments  and  proba¬ 
bility  distributions  of  acoustic  pressures,  and  how  these  are  affected  by  characteristics 
of  the  problem  such  as  wave  numbers. 

One  of  the  computational  tasks  required  is  the  solution  of  algebraic  systems  of 
equations  with  multiple  right  hand  sides.  In  the  case  of  uncertain  boundary  data,  the 
costs  of  this  component  of  the  computation  can  be  kept  low  using  the  fact  that  the 
solution  has  a  Kronecker  product  structure.  For  our  scattering  example,  the  systems 
can  be  solved  efficiently  with  a  multigrid  algorithm  for  the  discrete  Helmholtz  equation 
[9],  and  we  show  that  efficiency  can  be  enhanced  in  some  cases  using  block  iterative 
methods  for  systems  with  multiple  right  hand  sides  [4,  10,  17,  22].  With  this  strategy 
for  solving  the  algebraic  systems,  the  dominant  cost  of  the  computation  is  that  of 
computing  statistical  quantities.  We  also  show  that  the  Kronecker  product  structure 
of  the  solution  allows  storage  costs  to  be  kept  relatively  low,  and  moreover  it  enables 
the  statistical  computations  to  be  performed  using  efficient  matrix-oriented  operations 
that  are  trivially  parallelizable  and  amenable  to  implementation  using  Level  3  Basic 
Linear  Algebra  Subprograms  (BLAS3)  [8].  This  means  that  it  is  possible  to  handle 
relatively  fine  “discretization”  in  the  stochastic  domain  that  would  otherwise  not  be 
possible. 

We  note  that  an  alternative  approach  for  handling  random  right  hand  sides  has 
been  developed  in  Schwab  and  Todor  [20],  where  it  is  shown  that  the  mean  and 
second  moment  of  the  solution  can  be  computed  directly,  where  the  latter  entails  the 
solution  of  a  fourth  order  equation  derived  for  this  quantity.  It  is  shown  in  [20]  that 
when  the  underlying  differential  operator  is  coercive,  then  so  is  the  associated  fourth 
order  system,  and  efficient  multilevel  algorithms  (but  dependent  on  special  sparse 
grids)  can  be  developed  to  solve  it.  The  approach  under  consideration  here  has  the 
advantage  that  it  readily  yields  more  general  statistical  information  such  as  higher 
order  moments  and  probability  distributions.  It  is  also  relatively  straightforward  to 
implement,  essentially  only  requiring  algorithm  technology  for  second  order  problems. 
In  particular,  if,  as  in  the  example  considered  here,  the  underlying  problem  is  not 
coercive,  it  is  still  possible  to  take  advantage  of  efficient  algorithms  for  that  problem. 

A  summary  of  the  contents  of  the  paper  is  as  follows.  Section  2  contains  a 
description  of  the  stochastic  finite  element  methodology  and  identifies  the  structure 
of  the  algebraic  systems  derived  from  discretization.  Section  3  describes  the  iterative 
algorithms  that  we  consider  for  solving  the  discrete  Helmholtz  equation  and  the  block 
versions  designed  to  handle  multiple  right  hand  sides,  and  then  it  presents  some 
experimental  results  demonstrating  the  performance  of  these  solvers.  Section  4  then 
outlines  the  costs  of  computing  statistical  quantities  associated  with  the  solution  and 
shows  the  results  of  these  computations.  Finally,  Section  5  contains  some  concluding 
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remarks. 


2.  The  stochastic  finite  element  method.  We  briefly  describe  the  general 
methodology  with  an  eye  towards  showing  the  structure  of  the  algebraic  systems.  For 
our  description  we  use  the  problem  that  we  will  study  in  experiments,  the  Helmholtz 
equation;  it  will  be  obvious  that  the  approach  is  general.  See  [12]  for  complete  de¬ 
scriptions  of  this  methodology. 

2.1.  Introduction:  weak  formulation.  A  model  of  acoustic  scattering  from 
a  bounded  obstacle  is  given  by  the  Helmholtz  equation 

—A u  —  k2u  —  f  in  V 

B(u)  —  g  on  T  (2.1) 

t  =  L(u)  on  ^ 

where  the  solution  domain  V  C  Rd  is  bounded  internally  by  the  obstacle  boundary 

T  C  dV  and  externally  by  an  artificial  boundary  Y^.  The  boundary  differential 
operator  B  is  such  that  Dirichlet,  Neumann  or  Robin  boundary  conditions  result 
along  T,  and  L  is  the  Dirichlet-to-Neumann  operator  [13]  or  a  suitable  approximation 
thereof. 

The  weak  form  of  this  problem  is  to  find  u  E  Vg  such  that 

a(u,v)  =  £(v)  \/v  e  V  (2.2) 

where  V  and  Vg  denote  the  linear  and  affine  subspaces  of  H 1  (V)  of  functions  satisfying 
any  homogeneous  resp.  inhomogeneous  essential  boundary  conditions  along  T.  In 
the  simplest  case  of  Dirichlet  boundary  data  along  all  of  T,  the  sesquilinear  form 
a  :  H^V)  x  H^V)  ^  C  is 

a(u,v)  =  /  (Vu -Vv  —  k2uv)  dx  —  /  vL(u)  ds 

and  the  right  hand  side  functional  l  :  H 1  (V)  — >•  C  is 

£(v)  =  f  fv  dx. 

Jv 

To  introduce  randomness  into  this  formulation,  let  (II,  A,  P)  denote  a  probability 
space  with  sample  space  Yt ,  cr-algebra  A  and  probability  measure  P.  Let  (  :  II  — >•  C 
be  a  complex- valued  random  variable  with  £  E  L1(f I).  The  mean  or  expected  value  of 

C  is 


(0  =  [  CM  dP(uj)  =  [  z  dii{z), 

Jn  J  c 

where  g  is  the  distribution  probability  measure  associated  with  (  and  defined  on  the 
Borel  sets  B  in  the  complex  plane  by  g{B)  =  P((~1(B)).  Given  a  bounded  domain 
V  C  Rd  as  above,  a  random  function 

u  :V  x  i I  — >•  C,  (x,oj)  u{x ,  cj) 

is  one  that  is  jointly  measurable  with  respect  to  Lebesgue  measure  on  V  and  the 
measure  P  on  II  and  for  which 


(IK'MIMd))  <  oo. 
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The  space  of  random  functions  is  a  Hilbert  space  L2(V  x  Q)  with  respect  to  the  inner 
product 

{u,v)L2  =  {{u{x,  -),v(x,  ■))L2{V)). 

The  stochastic  Sobolev  spaces  Hk(V  x  ft)  are  defined  analogously. 

If  any  of  the  data  in  the  Helmholtz  equation  (2.1)  is  random  (e.g.,  the  wave 
number  k ,  forcing  function  /,  or  Dirichlet  boundary  data  g ),  then  the  solution  u  will 
be  a  random  function.  The  weak  form  of  the  stochastic  problem  is  then  to  find  u  E  Vg 
such  that 

(a(u,v))  =  (£(v))  Wv  e  Vo,  (2.3) 

where  Vg  and  Vo  are  the  stochastic  Sobolev  spaces  analogous  to  Vg  and  Vo- 

2.2.  The  Karhunen-Loeve  expansion  and  derived  weak  form.  We  con¬ 
sider  the  development  of  the  stochastic  finite  element  method  using  the  Karhunen- 
Loeve  (KL)  expansion ,  a  representation  of  random  functions  in  series  form  using  the 
eigenfunctions  of  the  covariance  operator.  For  the  sake  of  concreteness,  we  describe 
its  use  under  the  assumption  that  the  forcing  function  /  of  (2.1)  is  random;  we  will 
discuss  other  possibilities  in  Section  2.3. 

Let  the  covariance  function  associated  with  /  be  denoted  by 

c(x,y)  =  (f(x)f(y))  -  (f(x))  (f(y)) . 

Consider  the  integral  equation 

(C j)(x)  =  Xj(x),  where  (Cj)(x)  =  /  c(x,  y)j(y)  dy.  (2.4) 

Jt> 

This  is  a  linear  integral  eigenvalue  problem  in  which,  by  definition,  the  kernel  is 
symmetric  and  positive-semidefinite.  It  follows  from  the  general  theory  of  integral 
equations  [5,  Ch.  3]  that  C  is  a  compact  operator  and  there  exists  a  countable  se¬ 
quence  of  eigenpairs  {(A r,/r)}  where  the  eigenvalues  {Ar}  are  nonnegative  and  the 
eigenfunctions  {fr}  are  orthogonal  in  L2(V).  Let  the  eigenvalues  be  ordered  so  that 
Ai  >  A2  >  •  •  •  >  0.  The  Karhunen-Loeve  expansion  for  /  is 

00 

f(x,€)  =  fo(x)  +  \/Kfr(x)€r,  (2-5) 

r—1 

where  fo(x)  =  (/(#))  is  the  mean  of  /,  and  {£r(^)}r>i  are  uncorrelated  random 
variables  with  mean  zero  and  variance  one  [23,  pp.  447ff|. 

For  computation,  the  infinite  series  (2.5)  is  approximated  by  a  finite  one  with,  say, 
m  terms.  In  general,  the  more  localized  the  covariance  kernel  of  /  (the  smaller  the 
correlation  length),  the  slower  the  decay  of  its  eigenvalues  and  the  more  terms  need 
be  retained  in  the  KL  expansion  to  achieve  good  accuracy.  Thus,  the  utility  of  this 
approach  depends  on  the  assumption  that  the  properties  of  physical  systems  under 
consideration  vary  smoothly,  i.e. ,  there  are  significant  correlations  in  the  random  data. 
In  this  case,  it  is  expected  that  a  truncated  version  of  (2.5)  with  a  small  number  m 
of  terms  in  the  sum  is  sufficient  to  capture  the  randomness  in  the  system. 

Assume  now  that  the  random  function  is  given  by  a  finite-term  KL  expansion 

m 

f(x,0  =  fo(x)  +  ^2  (2-6) 

r— 1 
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Let  Xr  —  £r(fi)  denote  the  image  of  £r,  and  let  X  —  X\  x  •  •  •  x  Xm.  Collecting  these 
variables  into  the  random  vector  £  =  (£1? . . . ,  £m),  we  have  £(0.)  C  X.  Assume  that 
possesses  the  probability  density  function  pr  :  Xr  M,  which  gives  rise  to  the  joint 
density  function 


P{£)  =  Pl(6)P2(6)  •  •  •  Pm{im)- 

The  stochastic  variational  formulation  of  the  Helmholtz  equation  (2.1)  uses  as 
test  functions  random  functions  in  the  space 

V  =  ju(x,£)  :  (||u||/fi(x>))  p{£)d£,  <  ooj  ,  (2.7) 

with  trial  functions  in  the  space  Vg  defined  analogously.  The  stochastic  variational 
problem  is  then  specified  as  in  (2.3)  with 

{a{u-v)}  = 

{£{v))  =  (^  fv  dx  ^  p(£)d£. 

The  weak  solution  u  can  be  viewed  as  defined  on  a  (d  +  ra)-dimensional  domain  VxX. 

2.3.  Discretization  and  the  stochastic  system.  In  order  to  establish  no¬ 
tation,  we  briefly  discuss  the  discretization  of  the  deterministic  problem  (2.1),  as¬ 
suming  Dirichlet  boundary  conditions  u  =  g  hold  on  the  obstacle  boundary  T.  Let 
Vh  =  spanj^i, . . . ,  4>nx}  denote  a  finite  dimensional  subspace  of  Hq(T>),  and  let  Vg 
denote  the  affine  space  obtained  by  adding  basis  functions  . . . ,  (t>Nx+NE}  t° 

handle  degrees  of  freedom  on  the  boundary.  As  is  well  known,  the  discrete  weak 
formulation  entails  finding 


Nx  Nx+Ne 

Uh  =  Y,  U3(t>j  +  X  d(Xj)4>j 

3  — 1  j=Nx+ 1 


such  that 

Nx  r  Nx+Ne 

=  f<t>idx  -  Y  a(<l>j><l>i)9(xj)  Vi  =  l,...,Nx. 

i= i  j=Nx+ 1 

This  is  a  linear  system  of  equations  Au  =  f  where 

f  =  [(/,  -  AUEg;  (2.9) 

Ajje  represents  the  coupling  between  degrees  of  freedom  constrained  by  Dirichlet 
boundary  conditions  and  other  unknowns,  and  g  =  [g(xj)]fx^+ \  is  the  vector  of 
nodal  boundary  values. 

Now  consider  the  stochastic  problem  defined  by  (2.3)  and  (2.8).  For  the  dis¬ 
cretization,  let 

Vh  =  span{xj9(a;,£)  =  ^j{x)ipq{i)  :  j  =  1,. . .  ,NX,  q  =  1,...,^}, 
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denote  a  finite-dimensional  subspace  of  V  of  (2.7),  where  {?/q, . . . ,  }  is  a  basis 

for  a  finite-dimensional  subspace  of  L2(X).  Let  Vg  denote  the  affine  space  satisfying 
inhomogeneous  essential  boundary  conditions.  The  discrete  stochastic  problem  is  then 
to  find  Uh  £  Vg, 


N*  Nx 

uh(x,o  =  J2J2^j(x)ipq(0ujq 

q— i  j— l 


Nx+Ne 

+  ^2  fa(x)d(xj) 

j=Nx  + 1 


(2.10) 


such  that 


(a(uh,vh)}  =  (£(vh)}  V?;ft  G  Vh. 

The  result  is  a  linear  system  of  equations,  the  stochastic  system , 

Au=b  (2.11) 


of  order  Nx  x  A^,  for  unknowns 

(^11,  U2 1 ,  .  .  .  ,  Unx-1,N£,  UNx,N^)T  • 

Once  u  is  obtained,  statistical  properties  of  the  associated  random  function  Uh  can 
be  computed  easily,  see  Section  4. 

As  we  have  noted,  this  study  concerns  the  case  where  randomness  only  affects 
the  right  hand  side  of  the  algebraic  systems  generated,  i.e. ,  where  the  source  term 
or  boundary  data  is  random.  Let  us  consider  the  structure  of  the  discrete  problem 
(2.11)  in  this  case.  The  entries  of  the  finite  element  system  matrix  A  are 

(■ a(Xjq,XiP ))  =  f  P(€)d$ 

=  {^J  (yj  y<t>j -y4>i  -  k2(j)j4>idx  -  J  4>iL4>jdsj 

—  i'lpq'lpp)  5  $i)  5 

for  1  <  i,  j  <  Nx,  1  <  p,  q  <  N%.  Denoting  by  G  £  xNz  the  Grammian  matrix 

[G\Pq  =  (MP),  =  (2.12) 

and  by  A  £  CNx  xNx  the  stiffness  matrix  of  the  deterministic  equation,  the  coefficient 
matrix  is  seen  to  have  the  Kronecker  structure 

A  =  G®A. 

Note  that  this  implicitly  determines  an  ordering  for  the  rows  and  columns  of  A.  The 
rows  are  ordered  so  that  for  each  p,  indices  i  =  1, . . . ,  Nx  are  grouped  together,  and 
then  p  is  ordered  from  1  to  the  same  grouping  applies  to  the  columns. 

For  the  right  hand  side,  assume  as  in  Section  2.2  that  the  forcing  function  is 
random,  and  also  assume  for  the  moment  that  homogeneous  Dirichlet  boundary  con¬ 
ditions  g  —  0  hold  on  T.  It  then  follows  from  (2.6)  and  (2.8)  that  the  entry  of  f 
corresponding  to  the  test  function  \ip  —  ^i^p  is 

m 

(Z( Xip ))  =  /  ttf,  Xip)p{i)d£  =  t(fo,  <f>i)  (i’p)  +  L!  £(/r>  fa)  (tri^p)  •  (2.13) 

r—l 
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Let  us  define  the  vectors 


fr 

=  [(/rVi)]£l, 

r  =  0,1,.. 

. ,  m 

V’o 

=  [«V>]^1 

Vv 

II 

A 

r  =  0,1,.. 

. ,  m 

whereupon  the  discrete  system  has  the  form 

m 

(' G®A)u=f ,  f  =  -00  0  /q  +  yfc  (lpr  <g>  fr )  ■ 

r=l 


(2.14) 


That  is,  the  right  hand  side  lies  in  an  (m  +  l)-dimensional  subspace  of  xNx .  The 
solution  is  then 


u  =  (G®A)-1/  =  (G-1®T1)/ 

m  (2  15) 

=  (G- Vo)  ®  (^“7o)  +  ^  x/V  (G- vr)  ®  (A^fr) . 

r— 1 

This  entails  the  solution  of  m  +  1  systems  of  size  with  coefficient  matrix  G,  and 
m  +  1  systems  of  size  with  coefficient  matrix  A.  In  practice,  the  basis  functions 
{ipp}  for  the  stochastic  component  are  often  chosen  to  be  orthogonal  with  respect  to 
the  probability  measure  [7,  12],  in  which  case  G  is  a  diagonal  matrix.  Thus,  the  main 
computational  requirement  is  solution  of  the  m  + 1  systems  with  coefficient  matrix  A. 

Although  the  derivation  above  is  for  the  case  of  stochastic  forcing  function  and 
homogeneous  boundary  conditions,  the  conclusion  reached  is  general.  For  example,  if 
a  nonzero  Dirichlet  condition  holds  on  T,  then  the  construction  is  identical  except  /o 
has  the  form  (cf.  (2.9)) 

/o  =  [(/o,0i)]£i  -  AUEg- 

More  generally,  if  it  is  Dirichlet  boundary  conditions  that  are  random  (we  will  explore 
this  in  experiments  described  in  Section  3),  then  terms  of  the  form 

'ipQ  0  ( AUEgo )  +  ^2  Wv  ®  ( AEEgr )) 

r 

will  be  included  in  the  right  hand  side.  Similar  considerations  apply  for  Neumann 
conditions  on  the  obstacle  boundary. 

2.4.  Implementation.  The  notation  used  in  the  previous  section  treats  the 
unknowns  u  of  (2.11)  as  a  vector.  In  an  implementation,  it  is  in  fact  more  convenient 
to  treat  the  solution  as  a  two-dimensional  array.  In  particular,  consider  the  matrices 

F  [ fo  i  fl  i  i  i  fm\  5  A  dicig^  1 ,  \/  \\  ,  •  •  •  5  's/  ^  \^P  Qi^P  \i  5  •  •  •  i  ^P  m\i 

where  the  vectors  {fr}  and  are  defined  in  (2.14).  Then  the  system  (2.11)  is 

essentially  of  the  form 

AU  =  B ,  (2.16) 

where  B  =  FWT  with  W  =  G_1(4rA).  The  solution  can  then  be  represented  implic¬ 
itly  in  outer-product  form  as 

U  =  VWT,  (2.17) 

where  V  =  A~XF  is  obtained  by  solving  the  system  of  equations  AV  =  F  with  m  +  1 
right  hand  sides. 
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3.  Iterative  solution  of  the  stochastic  system.  For  the  problem  under  con¬ 
sideration,  the  coefficient  matrix  of  (2.16)  is  a  discrete  Helmholtz  operator,  which  is 
complex,  symmetric  and  indefinite.  In  this  section,  we  describe  an  iterative  algorithm 
that  can  be  used  to  solve  this  system  and  demonstrate  its  effectiveness  on  a  set  of 
benchmark  problems. 

3.1.  Solution  algorithm.  The  basic  solution  algorithm  we  use  is  a  multigrid 
method  designed  for  the  Helmholtz  equation,  adapted  to  handle  multiple  right  hand 
sides.  As  is  well  known,  the  principle  behind  multigrid  is  to  combine  smoothers  to 
eliminate  oscillatory  components  of  the  error  on  fine  grids,  together  with  coarse  grid 
corrections  to  eliminate  smooth  components.  For  the  Helmholtz  equation,  standard 
multigrid  approaches  are  not  effective.  There  are  two  difficulties: 

1.  Standard  smoothers  such  as  the  Jacobi  and  Gauss-Seidel  methods  do  not 
work  because  certain  smooth  modes  are  amplified  by  these  operations. 

2.  The  eigenvalues  associated  with  some  smooth  modes  change  signs  during 
the  grid  coarsening  process,  which  causes  the  coarse  grid  correction  to  also 
amplify  some  smooth  modes  rather  than  eliminate  them  from  the  error. 

These  difficulties  derive  from  the  indefiniteness  of  the  system.  In  [9],  we  developed 
a  method  that  addresses  them.  The  first  difficulty  is  handled  by  replacing  standard 
smoothers  with  Krylov  subspace  methods,  i.e. ,  GMRES  iteration  [19]  is  used  as  the 
smoother.  The  second  one  is  handled  by  using  the  multigrid  operation  as  a  precon¬ 
ditioner  for  an  outer  Krylov  subspace  iteration,  so  that  components  of  the  error  not 
treated  correctly  by  the  multigrid  coarse  grid  computations  are  eliminated.  Because 
the  multigrid  smoother  is  no  longer  a  linear  operator,  the  outer  iteration  must  han¬ 
dle  this  via  a  so-called  “flexible”  GMRES  algorithm  [18].  A  complete  description 
and  analysis  of  the  preconditioning  strategy  is  given  in  [9],  where  it  is  demonstrated 
that  the  algorithm  exhibits  “textbook  multigrid”  convergence  behavior,  that  is,  con¬ 
vergence  rates  that  are  independent  of  the  discretization  parameter;  there  is  some 
dependence  on  the  wave  number  k. 

We  also  adapt  this  approach  to  handle  the  system  (2.16)  with  multiple  right  hand 
sides,  the  number  of  which  is  denoted  by  m  within  this  section.  Recall  that  Krylov 
subspace  methods  generate  an  iterate  at  step  s  using  a  certain  subspace  of  dimension 
s.  Two  types  of  Krylov  subspace  algorithms  have  been  proposed  for  problems  with 
multiple  right  hand  sides: 

•  Block  algorithms  [2],  [17]  construct  a  subspace  of  dimension  ms  formed  by  the 
union  of  the  s-dimensional  subspaces  for  each  right  hand  side.  Then,  for  each 
right  hand  side,  they  find  the  best  solution  within  that  subspace.  Deflation 
is  used  to  remove  vectors  that  become  linearly  dependent. 

•  Seed  algorithms  [3],  [21]  form  a  Krylov  subspace  using  one  of  the  right  hand 
sides  and  then  find  the  best  solution  for  each  of  the  m  problems  within  that 
subspace.  If  the  seed  problem  converges  before  the  others,  then  a  different 
right  hand  side  is  chosen  as  the  seed  and  the  algorithm  is  repeated. 

Each  of  these  approaches  has  its  advantages.  Seed  methods  tend  to  perform  best 
when  the  right  hand  sides  are  related  to  each  other,  for  example,  if  they  arise  from 
functions  evaluated  at  nearby  points  [3].  This  approach  requires  less  storage  than 
block  methods:  for  systems  of  order  N ,  the  seed  GMRES  method  requires  storage 
proportional  to  sN ,  compared  to  smN  for  a  block  GMRES  solver.  On  the  other 
hand,  block  algorithms  tend  to  converge  more  rapidly  for  more  general  right  hand 
sides,  or  when  a  small  number  of  eigenvalues  are  well-separated  from  the  others  [17]. 
The  block  algorithm  also  makes  much  better  use  of  computer  memory  traffic,  since 
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each  access  to  the  coefficient  matrix  is  used  for  m  matrix- vector  products. 

In  our  application,  the  right  hand  side  vectors  (columns  of  F  in  (2.16))  derive  from 
the  orthogonal  eigenvectors  of  the  covariance  matrix,  and  we  found  the  seed  method 
to  be  ineffective.  Therefore,  we  restrict  our  attention  to  a  block  method.  The  idea  of 
block  iterative  methods  is  due  to  O’Leary  and  R.  Underwood.  The  block  biconjugate 
gradient  algorithm  was  described  in  [17],  and  a  block  quasi-minimum  residual  method 
in  [10].  Algorithms  for  altering  the  block  size  adaptively  were  given  in  [1].  A  block 
GMRES  algorithm  was  presented  by  Vital  [22]. 

We  also  need  to  modify  the  algorithm  to  handle  the  nonlinear  preconditioner, 
as  described  in  [4].  To  present  this  block  flexible  GMRES  method  for  (2.16),  we  use 
the  generic  notation  Ax  =  b  for  the  linear  system,  and  w  =  M(v)  to  represent 
a  generic  preconditioning  operation.  This  may  be  a  linear  operation  derived  from 
a  matrix,  or  (as  in  the  present  setting)  a  nonlinear  operation.  Let  denote  an 
approximation  of  the  solution  to  the  jth  equation  (of  m)  computed  at  iteration  s.  The 
block  flexible  GMRES  algorithm  generates  a  sequence  of  matrices  {Vj}  of  dimensions 
N  xm  that  together  form  a  matrix  V  =  [Vi , . . . ,  Vs] ,  and  a  set  of  matrices  Zj  =  M(Vj), 
and  Z  =  [Zi, . . . ,  Zs\.  The  block-Hessenberg  matrix  H  has  block  entries  Hij.  Each 
component  {x^  :  j  =  1, . . . ,  m}  of  the  iterative  solution  has  the  form 

x ^  =  xf^  +  Zy j8"1  (3.1) 

such  that  the  norm  of  the  residual  || bj  —  Ax^||  is  minimal,  where  x^  is  a  possibly 
arbitrary  initial  value.  If  M(v)  is  a  linear  operation  (say  Mv  where  M  is  a  matrix), 
then  Vs,  the  span  of  the  columns  of  V,  is  the  same  as 

span{Mru(MA)Mru. . . ,  (MAfl^Mn,. . . ,  Mrm,  (MA)Mrm, . . . ,  (M4)s_1Mrm}, 

the  span  of  the  Krylov  vectors.  In  this  case,  it  is  not  necessary  to  store  the  auxiliary 
matrix  Z.  If  M  is  a  nonlinear  operator,  then  Vs  will  not  be  a  Krylov  subspace,  but  the 
s’th  iterate  is  still  optimal  among  all  possibilities  of  the  form  (3.1),  which  corresponds 
to  an  affine  subspace  of  CN  of  dimension  ms.  The  algorithm  is  as  follows: 

Compute  the  residual  r  =  b  —  Ax  of  dimension  N  x  m. 
until  (\\ri\\  <  S,  £  =  1, . . .  ,ra), 

%  Generate  a  subspace  of  dimension  ms  from  the  residual  r. 

Define  V\  to  be  the  orthogonal  factor  in  the  QR  factorization  of  r. 
for  i  =  2, . . . ,  s  +  1, 

%  Generate  the  directions  defining  the  new  basis  vectors. 

Zi-i  =  M(U_!) 

W  = 

%  Orthogonalize  these  directions  against  the  previous  ones, 
for  j  =  l,...,i  — 1, 

Hpi-i  =  v;w 

TV  =  IV  -  VjHj'i- 1 
end  for  j. 

Perform  a  QR  factorization  of  TV,  obtaining  the  upper  tri¬ 
angular  factor  Hij-i  and  the  orthogonal  matrix  V\. 
end  for  i. 

%  Update  each  of  the  solutions, 
for  j  =  1  ,...,m, 
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Fig.  3.1.  Spatial  domain  and  initial  mesh  used  in  spatial  discretization. 
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c  =  V*rj 

Solve  the  least  squares  problem  min^  \\c  —  Hy\\ 
Xj  =  Xj  +  Zy 


end  for  j. 


r  =  b  —  Ax 
end  until 

The  loop  on  i  can  break  down  if  the  matrix  W  becomes  rank-deficient.  In  this  case, 
we  reduce  the  size  of  the  block  by  dropping  the  dependent  columns,  updating  the 
solutions,  and  continuing  with  the  residuals  that  have  not  converged. 

3.2.  Experimental  results.  We  tested  the  performance  of  the  block  flexi¬ 
ble  GMRES  algorithm  for  solving  the  stochastic  Helmholtz  equation  on  the  two- 
dimensional  domain  V  consisting  of  the  complement  within  the  unit  circle  of  a  scat- 
terer  taken  to  be  a  semi-open  cavity.  Dirichlet-to-Neumann  conditions  are  specified 
on  the  external  boundary  Too.  The  discretization  in  space  consists  of  piecewise  linear 
elements  on  triangles.  Figure  3.1  shows  the  scatterer  and  the  initial  mesh  used  in  all 
tests.  For  each  wave  number  k ,  this  mesh  is  refined  until  khmax  was  on  the  order  of 
7r/5  «  .63,  so  that  there  are  approximately  ten  points  per  wavelength.  All  compu¬ 
tations  were  done  using  Matlab.  Mesh  construction  was  done  using  the  Matlab 
PDE  TOOLBOX  routines  initmesh  and  ref  inemesh,  which  performs  a  uniform  mesh 
refinement. 

All  uncertainty  in  the  specification  of  the  boundary  value  problem  occurs  in  the 
statement  of  boundary  conditions  on  T,  the  boundary  of  the  scatterer,  where  Dirichlet 
boundary  conditions  u  =  g  are  such  that  g  is  a  random  function  as  specified  in  (2.6), 
with  mean  determined  by  an  incident  plane  wave  g(x  1,^2)  =  —elk^XlCOse+X2Sine^  at 
angle  0  =  7r/4.  We  assume  as  in  [6,  7]  that  the  random  variables  making  up  the 
KL  expansions  are  uniformly  distributed  on  an  interval  Zr  =  [—a,  a],  giving  rise  to 
the  joint  uniform  distribution  on  X  —  [—a,  a]m  with  joint  density  function 


(3.2) 


The  convention  that  (^j)  =  leads  to  the  condition  a  =  y/S.  The  stochastic 
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Table  3.1 

Number  of  matrix-vector  products  required  to  solve  m  +  1  systems  arising  from  m-term  KL- 
expansion,  using  preconditioned  FGMRES.  Numbers  in  parentheses  are  average  iteration  counts  or 
number  of  block  iterations. 


kh  =  .72,  Nx  =  4170 

m  —  4 

m  —  6 

m  =  8 

k  =  5tt 

Block  FGMRES 
FGMRES 

35  (7) 

37  (7.4) 

49  (7) 

52  (7.4) 

54  (6) 

67  (7.4) 

kh  =  .36,  Nx  =  16, 196 

m  =  4 

m  =  6 

m  =  8 

Block  FGMRES 
FGMRES 

40  (8) 

45  (9.0) 

56  (8) 

63  (9.0) 

72  (8) 

81  (9.0) 

kh  =  .72,  Nx  =  16, 196 

m  —  4 

m  —  6 

m  =  8 

k  =  107T 

Block  FGMRES 
FGMRES 

85  (17) 
153  (30.6) 

105  (15) 
214  (30.6) 

135  (15) 
276  (30.7) 

kh  =  .36,  Nx  =  63, 816 

m  —  4 

m  —  6 

m  =  8 

Block  FGMRES 
FGMRES 

90  (18) 
157  (31.4) 

119  (17) 
220(31.4) 

162  (18) 
282  (31.3) 

kh  =  .72,  Nx  =  63,816 

m  —  4 

m  —  6 

m  =  8 

k  =  20t r  Block  FGMRES 

FGMRES 

200  (40) 
360  (72.0) 

245  (35) 
495  (70.7) 

288  (32) 
636  (80.7) 

domain  X  —  [— is  discretized  using  a  uniform  mesh:  each  of  the  m  coordinate 
intervals  is  subdivided  into  n %  equal  subintervals,  resulting  in  =  n ™  elements,  each 
of  which  is  an  m-dimensional  cube  with  side  length  h %  =  2 a/n^.  Since  there  are  no 
continuity  requirements  on  the  probability  space,  the  basis  functions  are  taken  to  be 
piecewise  constants,  that  is,  the  basis  function  has  the  value  one  on  the  cube  with 
index  q  and  zero  elsewhere.  This  leads  to  the  particularly  simple  diagonal  structure 
for  the  Grammian  matrix  G  of  (2.12),  G  =  -^1.  (Recall  that  this  is  always  the  case 
if  the  basis  functions  for  the  stochastic  discretization  are  orthogonal  with  respect  to 
the  probability  measure.) 

To  define  the  KL  expansion,  the  covariance  function  associated  with  g  is  assumed 
to  have  the  form 


c(x !,x2)  =  a2e-( I(*i)i-(*2)iI)/2+I(*i)2-(*2)2I))>  xu  X2  g  r. 

A  general  requirement  of  this  methodology  is  that  the  first  m  eigenfunctions  and 
eigenvalues  of  the  covariance  operator,  or  discrete  approximations  to  them,  be  avail¬ 
able.  In  some  circumstances,  these  may  be  obtained  in  closed  form  [12,  pp.  27ff], 
or,  alternatively,  they  may  be  approximated  using  a  Galerkin  discretization  of  the 
integral  equation  (2.4).  For  the  problems  we  are  considering,  the  domain  of  g  is  one¬ 
dimensional,  and  the  Galerkin  computation  is  inexpensive.  The  discrete  eigenvalues 
and  eigenvectors  are  computed  directly  from  the  Galerkin  approximation.  If  the  do¬ 
main  in  question  is  of  higher  dimension,  this  computation  can  be  done  efficiently  using 
sparse  eigenvalue  methods  and  fast  summation  techniques  [14]. 
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Table  3.1  examines  the  performance  of  the  block  flexible  GMRES  algorithm  and 
compares  it  with  that  of  the  flexible  GMRES  algorithm  (FGMRES)  applied  to  each 
right  hand  side  separately.  In  these  tests,  the  stopping  criterion  for  the  solvers  was 
for  the  Euclidean  norm  of  each  component  of  the  residual  to  satisfy 

INI/IIM  <  10“6,  j  =  +  1. 

For  the  block  method,  the  iteration  was  stopped  when  the  maximal  individual  residual 
component  meets  this  criterion.  The  table  shows  the  total  number  of  matrix-vector 
products  performed  during  the  course  of  the  computation,  and  in  parentheses,  the 
number  of  iterations  required  for  convergence.  For  FGMRES,  the  latter  number  is 
the  average  for  m  +  1  right  hand  sides;  for  block  FGMRES,  it  is  the  number  of  block 
iterations.  Note  that  the  dimension  of  the  spaces  constructed  by  the  block  FGMRES 
method  depends  on  m,  the  number  of  right  hand  sides,  but  not  on  the  discretization 
parameter  n %  associated  with  the  stochastic  domain,  since  F  does  not  depend  on  n % 
in  (2.16). 

It  is  clear  from  these  results  that  the  block  methods  require  fewer  matrix- vector 
products  in  all  cases,  and  the  difference  in  the  number  of  these  products  becomes 
more  dramatic  as  the  number  of  right  hand  sides  increases  and  also  as  the  wave 
number  k  grows,  i.e. ,  as  the  problem  becomes  more  difficult.  The  results  provide 
further  evidence  of  the  mesh  independent  performance  of  the  multigrid  algorithm. 
We  note,  however,  that  as  the  number  of  steps  s  increases,  the  advantages  of  the 
block  FGMRES  method  become  less  pronounced,  since  the  overhead  in  generating 
the  Krylov  space  grow  like  m2s2NXj  compared  to  ms2Nx  when  the  right  hand  sides 
are  processed  separately.  A  block  Krylov  subspace  method  such  as  QMR  [10]  would 
not  suffer  from  this  drawback,  although  it  is  not  clear  that  this  approach  can  be 
adapted  to  handle  a  nonlinear  preconditioner.  Because  solution  of  the  linear  systems 
step  represents  a  low  order  cost  for  the  complete  construction  of  statistical  data  (see 
the  next  section),  we  have  not  explored  this  issue. 

4.  Computation  of  statistical  data.  Once  the  random  function  Uh  of  (2.10) 
is  available,  we  are  interested  in  statistical  properties  such  as  moments  and  probabil¬ 
ity  distributions  associated  with  it.  In  the  case  of  time-harmonic  wave  propagation, 
an  important  quantity  is  the  modulus  \uh |,  which  indicates  the  significance  of  the 
component  with  wave  number  k  in  the  wave  field.  In  this  section,  we  describe  the 
computations  required  to  generate  statistical  data  associated  with  the  random  func- 
tion  \uh(x,£,)\. 

4.1.  Computation  of  a  distribution  function.  Consider  the  construction  of 
the  probability  distribution  function  for  the  maximum  modulus 

F(a)  =  Pr(max  \v,h(x,£)\  <  a).  (4.1) 

X 


Let 


Sa  =  {£  el  :  max  |u/t(a;,£)|  <  a). 

X 

Using  the  definition  of  the  joint  density  function  (3.2),  we  have 

F(a)=Js  P(M=\Sa\(iy. 
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Fig.  4.1.  Effect  of  the  stochastic  discretization  parameter  n ^  on  the  estimate  of 
Pr( maxs  \uh(x,£)\  <  a),  for  k  =  and  various  choices  of  m,  the  number  of  KL  terms. 


•  -  --  — 1 - 1 - 1 - - —  o  *-»■■■«#* — _•  g  >  •-  4  i  -----  — j 
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To  determine  the  volume  of  Sa,  let  £  E  X  be  given,  and  let  q  =  q(£)  be  the  index  of 
the  stochastic  element  Tq  C  X  containing  £.  Then 


=  max  u 


'JQ  I 


max  \uh(x,  £)|  =  max 

X  X 

where  the  latter  equality  follows  from  the  linearity  of  Uh  in  space.  Letting 


Nx 


^2Ujq<Pj(x) 
3  =  1 


Sa 


{q  :  max  \  ujq\  <  a}  , 


it  follows  that  |<Sa|  =  sa,  and  therefore  F(a)  =  sa/N^.  This  construction  requires 
maxj  \ujq\  for  each  q.  Once  these  maxima  are  computed,  they  can  then  be  used  to 
compute  F(a)  for  any  a. 


13 


Fig.  4.2.  Effect  of  the  number  of  terms  m  in  the  KL  expansion  on  the  estimate  of 
Pr( maxs  \uh(xX)\  <  a)>  for  k  —  ^7r- 


We  show  in  Figures  4. 1-4.2  the  results  of  computing  the  distribution  function 
(4.1)  for  various  parameter  values.  There  is  no  analytic  expression  for  this  quantity, 
so  it  is  difficult  to  make  a  rigorous  assessment  of  the  accuracy  of  the  computations. 
Nevertheless,  it  is  possible  to  identify  certain  qualitative  aspects  of  the  results  as  well 
as  to  place  the  costs  of  producing  them  in  context.  First,  note  that  discretization  of 
the  random  component  of  the  problem  can  be  viewed  from  two  perspectives,  derived 
from  the  number  of  terms  m  used  in  the  finite  KL  expansion,  and,  once  m  is  fixed,  from 
the  value  n %  of  the  discretization  parameter  in  X.  Convergence  of  the  KL  expansion 
depends  on  the  correlations  within  the  process;  when  the  finite  expansion  is  fixed,  it 
is  shown  in  [7]  that  the  error  in  the  stochastic  discretization  (assuming  an  accurate 
spatial  discretization)  is  proportional  to  n^1 .  Since  the  number  of  stochastic  degrees 
of  freedom  is  proportional  to  =  n™,  it  would  be  desirable  to  keep  m  as  small  as 
possible. 

In  Figure  4.1,  we  consider  the  impact  of  the  two  parameters  m  and  n^,  for  a 
fixed  wave  number  k  =  57t.  (The  spatial  discretization  was  such  that  kh  =  .36.)  Each 
subplot  in  this  figure  corresponds  to  a  fixed  value  of  m  for  which  n %  is  allowed  to  vary. 
Each  plot  shows  convergence  to  a  fixed  curve  with  refinement  in  n^,  as  expected.  It  is 
also  noteworthy  that  as  m  is  increased,  the  quality  of  the  solution  for  fixed  n %  appears 
to  improve.  (For  example,  the  solution  for  =  4  is  closer  to  the  limiting  value  for 
each  successive  choice  m  —  2,  4,  6.)  This  indicates  that  the  constants  associated  with 
the  error  bounds  are  smaller  as  m  increases.  Figure  4.2  explores  the  impact  of  m  more 
closely.  For  this  example,  the  results  suggest  that  m  =  8  is  an  appropriate  limiting 
value  for  the  number  of  terms  in  the  KL  expansion.  With  n %  =  4,  this  yields  65,  536 
spatial  degrees  of  freedom.1  The  combination  of  smaller  values  of  m  together  with 
large  n %  (e.g.,  m  —  2  and  n %  =  10,  yielding  1024  stochastic  degrees  of  freedom)  is  able 
convey  the  qualitative  structure  of  the  distribution  at  significantly  smaller  cost. 


:We  also  remark  that  for  m  =  10,  n^  =  4  was  the  largest  discretization  parameter  we  could  use 
in  our  Matlab  environment.  This  led  to  =  1,048,576  stochastic  degrees  of  freedom. 
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Fig.  4.3.  Estimated  probability  distribution  function  Pr( maxs  \uh(x,£)\  <  a),  for  various 
values  of  the  wave  number  k. 


Finally,  Figure  4.3  shows  the  estimated  distribution  function  (4.1)  for  different 
values  of  k.  These  results  suggest  that  this  probability  distribution  function  does  not 
vary  dramatically  as  the  wave  number  increases. 

4.2.  Computation  of  higher  moments.  For  examples  of  other  statistical  data 
to  be  computed,  consider  the  moments  of  the  modulus  of  Uh .  Let  a^(x,  £)  =  \uh(x,  £)|, 
and  let 

a(h\x)  =  {ah{x,-)v) ,  v  =  1,2,3,... 
denote  the  moments  of  a^.  We  have 


*h\x) 


J  I Uh(x,0\vp(Od^ 

Xfz 

P=  1  J±P 


NZ  Nx 

Yf  (  ]  MO 

q= 1  \j=1 


p(m 


57  E 


Nx 

y^  ujp(t)j  (x) 

3  = 1 


s  p= i 


Nx 

y^  ujp4>j{x) 

3  = 1 


This  is  straightforward  to  evaluate  once  the  coefficients  are  available.  In  par¬ 

ticular,  the  nodal  values  are 


(n)  _  (v 

If  =al 


\*i) 


1 


y^i  Uip\u , 


p= 1 
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Fig.  4.4.  Estimated  mean  /ik?  standard  deviation  ratio  (Jhltthi  and  skewness  of  \uh\,  for 
k  =  5n. 
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giving  the  piecewise  linear  interpolant  of  the  z/th  moment, 


3  = 1 


Nx 


The  computations  required  for  central  moments  ^  —  a{h  }  j  y  are  identical  in  struc¬ 

ture. 

We  examine  some  of  these  quantities  in  Figures  4. 4-4. 6.  Four  things  are  shown: 
the  mean  standard  deviation  cr^,  ratio  of  standard  deviation  to  mean,  and  scaled 
third  central  moment  (the  coefficient  of  skewness  [16]) 


of  ah .  The  data  used  for  these  plots  come  from  the  parameter  choices  m  =  8  for 
the  truncated  KL  expansion,  stochastic  discretization  parameter  =  4  and  spatial 
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Standard  deviation  *  mean,  k=1  On 
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Skewness,  k=1  On 
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discretization  satisfying  kh  =  .36  for  both  k  =  hir  and  107r  and  kh  =  .72  for  k  =  207r. 
Within  each  figure,  the  means  and  standard  deviations  are  displayed  using  the  same 
scalings.  The  magnitudes  of  the  standard  deviations  largely  mirror  those  of  the  means, 
and  there  is  virtually  no  difference  in  the  relative  sizes  of  these  quantities  for  different 
wave  numbers.  This  indicates  that  size  of  the  wave  number  k  will  not  have  a  significant 
impact  on  the  confidence  that  can  be  attributed  to  computed  mean  solutions.  The 
depictions  of  skewness  indicate  that  near  the  corner  singularities,  the  distributions 
tend  to  be  more  skewed  toward  the  right  (positive  direction)  with  respect  to  the 
mean,  and  inside  the  cavity  they  are  skewed  more  toward  the  left;  this  may  be  of  use 
in  identifying  the  shape  of  scatterers. 

Note  that  all  these  computations  require  the  complete  set  of  values  {ujq  :  j  = 
1, . . . ,  Nx,  q  —  1, . . . ,  A^},  which  are  obtained  from  (2.17)  as 

m 

'U'jq  —  ^  ^  VjrWqr  . 
r— 0 
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Fig.  4.6.  Estimated  mean  standard  deviation  ratio  crh/fi^,  and  skewness  of  \uh\,  for 
k  =  207T. 
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Consequently,  the  cost  is  of  order  0(NxN^)  and  these  computations  represent  the 
dominant  expense  of  the  process.  The  storage  costs  are  also  of  this  magnitude  but 
can  be  reduced  to  order  m  ma x(Nx,N^)  by  taking  advantage  of  the  outer-product 
representation  (2.17)  and  recomputing  Ujq  whenever  it  is  needed.  The  tradeoff  here 
is  a  (small)  additional  computational  expense  of  magnitude  0(mNxN ^).  This  makes 
it  feasible  to  handle  large  values  of  m  or  n ^  that  storage  restrictions  would  otherwise 
prevent. 

5.  Concluding  remarks.  Our  aim  in  this  work  was  to  carefully  outline  the  com¬ 
putational  issues  associated  with  implementing  the  stochastic  finite  element  method 
and  processing  the  results  for  a  model  of  acoustic  scattering,  where  uncertainty  is 
restricted  to  boundary  data.  We  have  shown  that  a  representation  of  the  solution  in 
outer  product  form  leads  to  significant  savings  in  storage  and  also  enables  the  rel¬ 
atively  inexpensive  computation  of  the  random  solution.  The  dominant  cost  comes 
from  postprocessing  of  the  solution  to  compute  statistical  data,  although  the  outer 
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product  form  in  this  setting  reduces  the  storage  overhead  of  these  computations.  Fi¬ 
nally,  we  note  that  if  uncertainty  appears  in  the  differential  operator  instead  of  the 
right  hand  side,  then  the  outer  product  formulation  of  the  stochastic  system  is  not 
available,  and  this  problem  would  be  more  costly  to  solve. 

Acknowledgement:  The  authors  wish  to  thank  Ivo  Babuska  for  introducing 
them  to  this  topic  and  for  many  helpful  discussions. 
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